module func_drv_mod; 
use numeric_grid;
implicit none;
logical, save :: FUNC_DRV_TRACE_ENTRY = .TRUE.;
logical, save :: CACL_DRV_P9_TRACE = .TRUE.;

!from p24 of Bickley MathGaz 25 (1941)
! NOTE!!! Be careful, fortran mtrx is by-column
real(R14), dimension(5,9), save :: DRV_P9 = RESHAPE(& 
    (/-109584D0, 322560D0, -564480D0,  752640D0, -705600D0 &
    ,   -5040D0, -64224D0,  141120D0, -141120D0,  117600D0 &
    ,     720D0, -11520D0,  -38304D0,   80640D0,  -50400D0 &
    ,    -240D0,   2880D0,  -20160D0,  -18144D0,   50400D0 &
    ,     144D0,  -1536D0,    8064D0,  -32256D0,       0D0 &
    ,    -144D0,   1440D0,   -6720D0,   20160D0,  -50400D0 &
    ,     240D0,  -2304D0,   10080D0,  -26880D0,   50400D0 &
    ,    -720D0,   6720D0,  -28224D0,   70560D0, -117600D0 &
    ,    5040D0, -46080D0,  188160D0, -451584D0,  705600D0/)&
    , (/5,9/));

contains 

function new_func_drv(px, py) result(res);    implicit none;  
  real(R14), dimension(:), pointer :: px, py, res;      
  call dbg_info("-->new_func_drv()", FUNC_DRV_TRACE_ENTRY);
  res => new_vec_from(py);
  call func_drv_p9(px, py, res);
  return;
end function

! 9-point derivative
! NOTE: this func fast but could be made faster if needed.
! Not tested for dimension(J:J2), where J != 1
SUBROUTINE func_drv_p9(px, py, res);    implicit none;  
  real(R14), dimension(:), pointer :: px, py, res;      
  real(R14) stx, h2; ! step x
  integer xn, mini, maxi, i, k, fromi;
  call dbg_info("-->func_drv_p9()", FUNC_DRV_TRACE_ENTRY);
  xn = SIZE(px, 1);  mini = LBOUND(px, 1); maxi = UBOUND(px, 1);
  call assert_true("xn >= 9", xn >= 9);
  call assert_true("xn == SIZE(py, 1)", xn == SIZE(py, 1));
  call assert_true("xn == SIZE(res, 1)", xn == SIZE(res, 1));
  call assert_true("mini == LBOUND(py, 1)", mini == LBOUND(py, 1));
  call assert_true("mini == LBOUND(res, 1)", mini == LBOUND(res, 1));
  call assert_step_grid_4(px);
  stx = px(mini+1) - px(mini);
  h2 = D05 / (20160D0 * stx);
 
  call dbg_vec("px=", px);
  call dbg_vec("py=", py);
  fromi = mini;
  do k = 1, 4 ! first 4
    i = fromi + k - 1;
    res(i) = h2 * calc_drv_p9(k, fromi, py);
    !call dbg_int("i=", i);
    !call dbg_real("pdy(i)=", pdy(i));
  end do
  k = 5;
  do i = mini+4, maxi-4 ! middle
    fromi = i - k + 1;    ! i = fromi + k - 1;
    res(i) = h2 * calc_drv_p9(k, fromi, py);
    !call dbg_int("i=", i);
    !call dbg_real("pdy(i)=", pdy(i));
  end do
  i = maxi-8;  
  do k = 6, 9 ! last 4
    i = fromi + k - 1;
    res(i) = h2 * calc_drv_p9(k, fromi, py);
    !call dbg_int("i=", (i+k-1));
    !call dbg_real("pdy(i)=", pdy(i));
  end do
  
  return;
end SUBROUTINE

real(R14) function calc_drv_p9(k, fromi, py) result(res);    implicit none;  
  real(R14), dimension(:), pointer :: py;      
  integer midi, k, fromi, mk;  
  logical :: trace;
  trace = CACL_DRV_P9_TRACE;
  !call dbg_int("k=", k, trace);
  !call dbg_int("fromi=", fromi, trace);
  midi = fromi + 4;   !call dbg_int("midi=", midi, trace);
  mk = 10 - k; !mirror k

  !call dbg_real("DRV_P9(5, k)=", DRV_P9(5, k), trace);  
  res = DRV_P9(5, k) * py(midi); ! middle
  !call dbg_real("res=", res, trace);  

  !call dbg_real("  DRV_P9(4, k)=", DRV_P9(4, k), trace);  
  !call dbg_real("DRV_P9(4, mk)=", DRV_P9(4, mk), trace);  
  res = res + DRV_P9(4,k)*py(midi-1) - DRV_P9(4,mk)*py(midi+1);
  !call dbg_real("res=", res, trace);  

  !call dbg_real("  DRV_P9(3, k)=", DRV_P9(3, k), trace);  
  !call dbg_real("DRV_P9(3, mk)=", DRV_P9(3, mk), trace);  
  res = res + DRV_P9(3,k)*py(midi-2) - DRV_P9(3,mk)*py(midi+2);
  !call dbg_real("res=", res, trace);  

  !call dbg_real("  DRV_P9(2, k)=", DRV_P9(2, k), trace);  
  !call dbg_real("DRV_P9(2, mk)=", DRV_P9(2, mk), trace);  
  res = res + DRV_P9(2,k)*py(midi-3) - DRV_P9(2,mk)*py(midi+3);
  !call dbg_real("res=", res, trace);  

  res = res + DRV_P9(1,k)*py(midi-4) - DRV_P9(1,mk)*py(midi+4);
  !call dbg_real("res=", res, trace);  
  return;
end function

end module 